Volumetric medical imaging techniques that can produce three-dimensional medical imaging data using any of a variety of imaging modalities, for example CT, PET, MRI, ultrasound, and X-ray, are now widely used for imaging or diagnostic purposes.
Volumetric medical imaging data may comprise a three-dimensional array of voxels, each voxel representative of a particular position in three-dimensional space and each voxel having one or more data values. For example in the case of CT data, each voxel may comprise an intensity value representative of the attenuation of the applied X-ray radiation provided at the location represented by the voxel, for example the attenuation of the part of a patient's anatomy present at that location.
A wide variety of rendering techniques are known that are used to process volumetric image data in order to display a desired representation of a patient or other subject. For example, the volumetric image data may be rendered to display a two-dimensional projection of the three dimensional volumetric data. The two-dimensional projection may include desired texture or shading effects, and may provide an impression of a surface in three dimensions to the viewer.
Commonly used volume rendering techniques use a transfer function to map each of a series of sampled voxel values to an opacity value and a colour value (usually represented by a combination of red, green and blue colour values). In ray casting and similar techniques, for each pixel in the 2D projection an imaginary ray is cast from the perspective of the viewer through a corresponding point on a two-dimensional imaginary image plane and through the volume represented by the volumetric image data.
The full volume rendering integral can be complex to evaluate accurately in real time. Therefore, it is known to fall back to a simple numerical approach of sampling c(s) (the colour of the sample) in discrete steps along rays cast through the volume, as illustrated schematically in FIG. 1. For a series of (usually equally spaced) sample points along the ray path through the volume, the transfer function can be used to assign opacity and colour values to the sampled points. A process is then performed to combine the colour and opacity values of the sampled points in sequence along the ray path to produce a single colour value for that pixel.
If, for example, the opacity of the first sampled point is high, as determined by the transfer function based on the measured intensity value for that point, then the colour of the corresponding pixel should be determined entirely, or almost entirely, by the colour value for that sampled point, as the imaginary ray effectively does not penetrate into the volume beyond that initial high opacity point. If the opacity values of the first points in the series of sampled points are relatively low then the colour of the corresponding pixel should be determined from a combination of the colour values for those sampled points in the series up to the point where the accumulated opacity reaches a sufficiently high value.
The simple sampling approach illustrated in FIG. 1 can be vulnerable to volatility in the transfer function and in the data sampled. An alternative known approach uses a pre-integration approach, which can go some way to account for volatility in the transfer function and in the data samples the renderer receives. Pre-integration is a technique that involves using colours that are pre-computed for the sample interval, with the results of the computations being stored in a 2-D look-up table (front-sample, back-sample) as illustrated schematically in FIGS. 2a and 2b. Corresponding opacity and colour values for each possible pair of intensity values and sample spacings are pre-computed using the transfer function and stored in the look-up table.
For each successive pair of samples along the ray path, indicated as front slice (Sf) and back slice (Sb) in FIGS. 2a and 2b, the renderer looks up the pre-computed colour and opacity value from the 2-D look-up table that corresponds to that pair of sampled intensity values and that sample spacing.
FIG. 2a is a representation of the transfer function, including red 10, green 12, blue 14 and opacity (α) 16 channels. The intensities of the particular sampled points, Sf and Sb, are indicated in FIG. 1a by dashed lines. Rather than calculating colours and opacities for each of the Sf and Sb separately and then including in a resulting integration process, in the pre-integration approach the opacity value and colour value for the pair of values at the sampled points Sf and Sb can be read directly from the 2D look-up table. The 2D look-up table is represented graphically in FIG. 2b as a 2D intensity plot. Effectively, the approach presumes that for each pair of sampled intensity values, the intensity varies linearly (or in some other specified way) in the interval between the sampled points, and thus the stored value in the look-up table for a particular pair of sampled intensity values can he considered to represent averaged opacity and colour values of the transfer function for the range of intensities between the two sampled intensity values.
The use of the 2D pre-integration can turn the estimate provided by the sampling approach of FIG. 1 (shown again in FIG. 3a) to a smoother estimate, as illustrated schematically in FIG. 3b, although very noisy data can still cause inaccuracies. Effectively, each block is now replaced by a pre-computed integral of the transfer function (between the two samples). This computation is subject only to the sampling rate and the transfer function, independent of the volume data, the viewing angle and all other rendering parameters. The 2D pre-integration approach effectively provides a slab-by-slab rendering process rather than slice-by-slice rendering process for each pixel, and thus can provide a smoothing to counteract volatility in the transfer function and in the data sampled.
The simple numerical sampling approach, for example as described with reference to FIG. 1, can lead to artefacts such as the tree ring artefacts shown in FIG. 4a that results from a thin skin layer going un-sampled in particular areas due to the discretization of the simple numerical approach. The 2D pre-integration approach captures the skin detail more accurately, as shown in FIG. 4b. 
One limitation of the 2D pre-integration approach is that precise 2D pre-integration can consume too much memory. For example, for the range of a 16 bit unsigned integer a floating point RGBA representation costs 128 bits×655362=64 gigabytes. Even before combining that requirement with the memory and processing requirements associated with any segmentation (which may require one transfer function per object) that usually represents an unacceptable burden.
The compute time for 2D pre-integration approaches can also be long. The nature of the integral involves many slow-to-compute procedures such as power functions (xy). This computation also has to be redone for every change in the sampling rate, usually left to the control of the application developer.
For many practical purposes, such as medical imaging applications for use by medical practitioners, speed and memory requirements are significant factors. Medical practitioners, such as doctors, radiologists, nurses and radiographers, do not expect to wait significant periods of time for rendered images to be produced, indeed they may expect such renderings to be produced almost instantaneously using processing resources available on a desktop terminal. Furthermore, medical practitioners commonly wish to change rendering parameters, leading to changes in transfer functions, or navigate through images whilst they are viewing. Again, the medical practitioner does not usually expect to wait for any significant periods of time for new renderings to be produced following such changes.
It has been proposed to adopt a 1D pre-integration approach, based on a 1D look-up table to approximate the functionality of the 2D pre-integration approach whilst reducing memory and processing burdens. One example of such a 1D approach is described in H. Kye et al, Interactive Classification for Pre-Integrated Volume Rendering of High-Precision Volume Data, Graphical Models, 70 (2008), 125-132.
In 1D approaches such as that of H. Kye et al, the look-up table contains values of colour and opacity that are pre-computed as corresponding to each possible value of measured image data (for example, each possible measured intensity value). In some cases the 1D table may have lower resolution than the measured image data (for example it may store entries only for integer intensity values) but stored values may be interpolated if desired. Effectively, the 1D table may be considered to store colour and opacity values pre-computed using the transfer function for pairs of sampled intensity values, similar to the 2D approach, but in this case one of the intensity values of the pair is set always to be at zero or some other common value.
It should be noted that in some cases the 1D look-up table may consist of two associated 1D look-up tables, for example two interleaved 1D look-up tables, which one table stores a series of cumulative integrals between 0 and a series of image data values (e.g. intensity values) and the other table stores the reverse cumulative integrals between infinity (or a maximum value) and the series of image data values (e.g. intensity values).
The 1D approaches described in the preceding two paragraphs can fit within memory limitations as it is only necessary to store a table as large as the classifier being used. Computation times are also generally shorter than with the 2D approach as the computations involved are summations rather than exponentials, which are faster to compute. The image quality provided by the 1D approaches can also be better than that provided by the simple sampling approach of FIG. 1.
However, known 1d approaches such as that of H. Kye et al generally use an approximation of the opacity of the sampled colour which is inaccurate for high opacity. The graphs of FIGS. 5a and 5b are contour plots of estimated opacity as a function of original opacity (x-axis) and ray step-size (y-axis) for accurate step-size weighted opacity (FIG. 5a) and for the approximation used in the 1D method (FIG. 5b). It can be seen that for values of opacity close to 1.0 (approaching complete opacity) the plots are more different, i.e. the approximation used in the 1D method is worse and providing too high transparency (too low opacity) compared to the true situation.
The effects of the faulty opacity approximation for high opacity values can be seen, for example, in FIG. 6 of H. Kye et al, Interactive Classification for Pre-Integrated Volume Rendering of High-Precision Volume Data, Graphical Models, 70 (2008), 125-132, in which the transfer function is set up to render the skin of a subject and the 1D method renders the skin more transparent than it should be, resulting in a different view that, undesirably in this case, renders internal structure such as bone. The computed opacity for a region in 1D methods may, for example, be representative of an average of the opacity within the region—thus a region containing an opaque interval and a non-opaque interval is assessed as non-opaque.
The faulty opacity approximation for high opacity values in 1D methods can also, more subtly, result in a different colour being rendered than should be the case if the opacity was approximated more accurately. That effect can be understood with reference to FIG. 6 of the present application, which represents a transfer function as a plot of opacity 20 and red 22, green 24 and blue 26 colour values as a function of intensity. In this case, if at one of the sample points (indicated by arrow 30 in FIG. 6) the value of intensity from the image data had an intensity value slightly greater than the value of intensity at which opacity first reaches a value of 1.0 according to the transfer function (indicated by arrow 28 in FIG. 6) then the resulting colour produced will have greater green and blue values than should be the case.
As mentioned above, in the 1D approach the pre-computed values of a opacity, and red, green and blue stored in the look up table for a particular sampled value of intensity are effectively an average of the values of opacity and red, green and blue obtained from the transfer function for all intensity values up to the sampled value. There may also be stored in the look-up table pre-computed values of opacity, and red, green and blue stored that are effectively an average of those values from the transfer function for all intensity values from infinity or a maximum value down to the sampled value. In the example of FIG. 6, the stored opacity and red, green and blue values in the look up table for the intensity 30 of the sample point will include contributions from points on the transfer function after point 28 when the opacity first reaches its maximum value of 1.0. However, in reality light would not penetrate beyond point 28, as the subject is considered to be completely opaque at that point. Therefore, in the example of FIG. 6 the resulting colour will be more green and blue than expected due to the fact that the pre-rendering process effectively samples behind something that is opaque and continues to pick up erroneous colour contributions as a result.
An alternative approach uses adaptive sampling, in which sampling is performed at a higher rate in high gradient areas. However, this approach generally fails to capture high frequencies in the transfer function, and can also reduce performance much more than expected because it interferes with parallelisation/prediction models and other assumptions involved in current rendering techniques.
Known rendering processes often provide shading and lighting effects in the rendered images as well as setting colour or greyscale to an appropriate value for each pixel. Such shading and lighting effects are usually provided by specifying the location of one or more imaginary light sources and then, as part of the rendering process, determining the effect of the light source(s) on the brightness of each pixel in the rendered image. The appropriate shading, represented for example by pixel brightness level, for each pixel is usually determined by analysing surface orientation around each of the sample points along the ray path through the volume for that pixel. A gradient volume (for instance a 4×4×4 voxel volume) may be defined around each sample point and a known gradient determination process is performed based on the intensity values for each of the voxels in the gradient volume. The resulting gradient values for the sampled points are then included in the subsequent integral process that combines the colour, opacity and, in this case, gradient values of the sampled points in sequence along the ray path to produce a single colour and shading/brightness value for that pixel.
As noted above, known 1D pre-integration approaches to sampling can produce approximation errors that can cause the opacity approximations for high opacity values to be underestimated. Effectively, such opacity approximation errors can allow the imaginary ray to progress too far along the ray path, which can cause interference from near lying surfaces when determining local gradients for use in the shading/brightness rendering. Even small traversal errors can result in catastrophic shading errors in some cases.
Shading errors resulting from underestimation of opacity in 1D methods are illustrated by way of example in FIGS. 7a to 7c, which show images of a skull rendered using a classification method such as that described in relation to FIG. 1 (FIG. 7a), rendered using a 1D look-up table method (FIG. 7b) and rendered using a 2D look-up table method (FIG. 7c). Hand-drawn arrows in FIG. 7b point towards regions of the image that include artefacts comprising unwanted texture features due to the under-estimation of opacity in the 1D look-up table method. Magnified parts of the images are shown in FIGS. 8a to 8c. 
The errors in the rendered images due to the underestimation of opacity for high opacities in the 1D look-up table methods can also be understood with reference to FIGS. 9a to 9d. FIG. 9a is a schematic illustration showing a sample point 40a along the path 42 of an imaginary ray cast through a volume represented by imaging data, for example a set of voxels of CT imaging data. The imaging data comprises a region 44 of high-intensity (and high opacity according to the transfer function in this case) with all other regions being of low intensity (and thus low opacity according to the transfer function in this case). A gradient volume 46a around the sampling point 40a is indicated in FIG. 9a. Voxels within the gradient volume 46a are included in a gradient calculation process to determine the orientation of any surfaces. The gradient data resulting from the gradient calculation process is associated with the sample point 40a and is included in a subsequent integral process that combines the colour, opacity and, in this case, gradient values of a series of sampled points 40a-40d in sequence along the ray path 42 to produce a single colour and shading/brightness value for a pixel of the ultimate rendered image.
FIGS. 9b to 9d show a series of further sample points 40b-40d as the sampling moves along the ray path 42. For each of the sampling positions, opacity and colour data is obtained from the 1D look-up table, and gradient data is obtained by processing intensity data for voxels within the gradient volume 46b-46d for that sampling position.
It can be seen that the sample point 40c has reached the high opacity region 44 in FIG. 9c. The ray should, ideally, terminate at that point as, due to the high opacity of region 44, the accumulated opacity for the series of sample points is such that the ray should not be able to penetrate further (in this example, it can be noted that the accumulated opacity is equal to the opacity at point 40c due to the low opacity at the preceding points 40a, 40b).
However, as the 1D look-up table underestimates the opacity for high opacity regions, the sample point progresses to the next position 40d, as shown in FIG. 9d. At that sample point 40d a further iso-surface of the region 44 enters the gradient volume 46d, thus affecting the calculation of the gradient data for that sample position 40d and so causing a shading error in the pixel in the ultimate rendered image corresponding to that ray path 42.